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We study the linear phenomenological Maxwell's equations in the presence of a polarizable and 

magnetizable medium (magnetodielectric). For a dispersive, non- absorptive, medium with equal 

electric and magnetic permeabilities e(w) and fJ>(u>), the latter can assume the value —1 (+1 is their 

vacuum value) for a discrete set of frequencies ±a>„, i.e., for these frequencies the medium behaves 

as a negative index material (NIM). We show that such systems have a well-defined time evolution. 

In particular the fields remain square integrable (and the electromagnetic energy finite) if this is 

the case at some initial time. Next we turn to the Green's function G(x, y, z) (a tensor), associated 

with the electric Helmholtz operator, for a set of parallel layers filled with a material. We express it 

in terms of the well-known scalar s and p ones. For a half space filled with the material and with a 

single dispersive Lorentz form for e(uu) = fJ>(u>) we obtain an explicit form for G. We find the usual 

behavior for negative index materials for uj = ±w, there is no refection outside the evanescent regime 

and the transmission (refraction) shows the usual NIM behavior. We find that G has poles in ±d>, 

^ ' ' which lead to a modulation of the radiative decay probability of an excited atom. The formalism is 

^) , free from ambiguities in the sign of the refractive index. 

■I— I , 
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P PACS numbers: 03.50.De,78.20.Ci, 42.25.Gy 
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O 1 I. INTRODUCTION 

■ I— I " 

Often magnetization plays a minor role in situations where the phenomenological Maxwell's equations apply. But 
in recent years negative index materials (NIM's), also called left handed materials, have become of increasing interest, 
in particular due to the work of Veselago |l| and Pendry 0] ■ Here the magnetization is not negligible at all. In general 
a NIM system is defined by the property that for certain frequencies w the electric permeability (permittivity) e{oj) 
or the magnetic permeability /i(w) becomes negative. Of particular interest is the case where both become negative 
at the same frequency w, the NIM frequency, and are equal to the opposite of their vacuum value, i.e., —1 instead of 
+1. Below we refer to this case as the NIM situation. 

For an introduction, containing an extensive set of references, see [3| (there are also NIM systems based upon 
specific properties of photonic crystals, which are not considered here). 

The existence of NIM's has been debated in the theoretical literature at various occasions Q. In particular the sign 
of the index of refraction, which involves taking a square root, has been a subject of discussion. Naively it equals +1, 
in both vacuum and a NIM system but this result is challenged for the NIM situation. For experimental verification, 
see [H. Calculations based upon a simple model, where one part of space is vacuum [e = fi = 1) and the other filled 
with a NIM (e = [i = — 1, frequency- independent) , tend to give ambiguous results. This is sometimes remedied by 
adding a small imaginary part to one of the permeabilities but on the whole the situation is rather unclear. 

The use of the phenomenological Maxwell's equations should solve possible ambiguities but it seems that so far 
this approach has not been taken and here we intend to fill this gap. Since fabricated materials, intended to study 
NIM behavior, are usually anisotropic, we take the (space and frequency dependent) susceptibilities (which relate the 
polarization and magnetization to the electric and magnetic fields) to be tensors rather than scalars. 

The first matter to be solved is the existence of a proper time evolution. In view of the time convolutions in the 
constitutive equations this is not directly obvious. The next task is to see if a NIM situation can exist. This being 
the case, the following point of interest is obtaining the Helmholtz Green's function and scattering amplitudes for 
specific configurations. The former is important since its imaginary part enters the radiative decay rate of an atom or 
nanostructure close to the material. Hence experimental results on such decay rates can give information about the 
properties of the material. In addition the Green's function, or rather the associated transition operator, is required 
to describe scattering phenomena, such as reflection and transmission in layered systems. 

Thus we start with the phenomenological Maxwell's equations with general frequency-dependent permeability ten- 
sors, satisfying the usual causality and passivity conditions. After providing some relevant background and a summary 
of the properties of the electric and magnetic susceptibilities we introduce the auxiliary field formalism (AFF). The 
latter was presented earlier by one of us (AT in Q) for dielectrics (^ = 1). The idea is to introduce an additional set of 
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fields, the auxiliary fields, to remove the time convolutions in Maxwell's equations. This has a number of advantages: 

1) The combined set of electromagnetic and auxiliary fields satisfies a unitary time evolution, thus insuring a proper 
time evolution for the electromagnetic fields. 

2) Such a system is easily quantized, leading to a second quantization formalism that is rigorously valid for both 
absorptive and dispersive systems. 

3) The formalism implies that the inverses of the electric and magnetic Helmholtz operators exist as bounded operators, 
so the associated Green's functions are square integrable. 

4) Setting up a scattering formalism is straightforward. 

The AFF leads to a proper time evolution, notwithstanding the possibility that for specific frequencies u> we can 
have a NIM situation, e(d>) = n(ui) = — 1. In case the initial fields are square integrable they remain so for all later 
times. In the Appendix we give a rigorous proof of this important fundamental property. 

Another relevant piece of information is that the susceptibilities for general dispersive, non-absorptive, systems 
consist of a (possibly infinite) sum of Lorentz contributions (AT in Q). This immediately gives a positive answer to 
the existence of NIM's. In case we are dealing with a single dispersive Lorentz contribution 

e (o,)=^( w ) = l--A- j, (1.1) 

CO — UJq 

we note that for the NIM frequences u> — ±d), a) 2 = w 2 + f2 2 /2 we have e(±w) = /i(±cD) = — 1. Thus, theoretically, the 
NIM case can be realized for dispersive, non-absorptive, systems, contrary to what is sometimes claimed [![. Adding 
more Lorentz terms gives more frequency values with this property but it remains a discrete set and in between the 
values of e{ui) and vary wildly. In fact, between two subsequent NIM frequencies there is always a frequency for 
which total reflection takes place (for lu = uiq in the above example), precisely the opposite of the NIM case, where 
no reflection is thought to be the situation. Thus it seems that obtaining an extended frequency interval for which 
the permeabilities are approximately equal to —1 is not possible. In fact it was already noted by Veselago [l[ that a 
system showing NIM behavior must be dispersive. In case there is absorption, s(uj) = /i(w) = —1 cannot be realized 
for real to as can be seen by adding absorption to the above case, 

e(u)=n{u) = l- " (1-2) 

W z + IJUJ — U!q 

Next we introduce the Laplace-transformed Maxwell's equations and the tensor Green's function G(x, y, z) related to 
the electric Helmholtz operator. G(x,y, z) features most of the properties of the system. As said the radiative decay 
rate of excited atoms is proportional to its imaginary part. We then turn to layered systems and express G(x, y,z) 
into a set of two scalar ones, for s and p polarization, respectively. Subsequently the Green's function for the half 
space case, mentioned above, is studied. In particular we obtain an explicit expression for G(x, y, z). Then E(x,t) is 
given by the inverse Laplace transform of 

£(x,z) = / dyG(x, y, z) ■ g(y, z), (1.3) 



where g(y, z) is some square integrable initial field configuration or an external current density. At this point the 
square root of z 2 e(x, z)/i(x, z) — k 2 , with k a two-dimensional wave-vector, must be evaluated as S > in z = lo + iS 
tends to 0. Depending on the values of x and u> different results are obtained, it can be positive, negative or imaginary. 
We find that reflective contributions to G(x, y,±cD + iO) vanish in the radiative regime, u) > k, and transmission is 
also modified substantially. This confirms the results by Pendry Q. In the evanescent regime G(x, y,z) has poles 
in z — ±cD, giving finite contributions, proportional to exp[±id)i], to the electric field 22(x, i). It turns out that K, 
the generator of the time evolution in the AFF, has ±ui as eigenvalues with infinite degeneracy, the latter giving rise 
to the above poles. Although we do not discuss quantization, we note that this feature gives rise to an interesting 
structure of the associated field Hamiltonian. In addition to the eigenvalue 0, associated with the vacuum state, now 
±<D are also eigenvalues. This affects radiative decay constants of excited atoms, as is discussed in Section VII. 

A word about notation: With a dispersive system we mean a dispersive, non- absorptive system. Inner products are 
denoted as (f,g) — (g\f)- The unit vector along a € R 3 is e a = a/a, a = |a|. The three Cartesian axes are denoted 
by Xi, X-i and A3 with corresponding unit vectors ei, e.^ and 63. The component of a _L e3 is denoted by a^. U 
is the unit 3x3 matrix. Transposes of matrices are indicated by means of the superscript T and their Hermitean 
adjoints by f . Square roots are defined in the usual way with non- negative imaginary part. I^(x)is the characteristic 
function for the set A, \a(x) — 1 for x £ A and \a( x ) = for x ^ A. 
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II. BACKGROUND 



Starting point is the set of linear phenomenological Maxwell's equations for the case that permanent polarization 
and magnetization are absent (we set e = Mo = 1 f° r brevity) 

d t D{x,t) = 9 X x H(x,t), dtB(x,t) = -d x x B(x,t), 
ftfD(x,t)=0, fl K .B(x,t)=0 J (2.1) 

with the constitutive equations 

D(x, t) = E(x, t) + P(x, t), P(x, i) = / ds Xe (x, * - s) • £(x, s), 

■/to 

(x, t) = B(x, t) - M(x, t), M(x, i) = f ds Xm (x, t-s)- (x, s). (2.2) 

Jto 

Here X e ( x , £) and Xm( x , i) ar e the electric and magnetic susceptibility tensors. We also introduce the current densities 
J e (x, t) — 9 t P(x, t) and J m (x, t) = d t M(x, t). Causality requires that the susceptibilities X e ( x , t) and Xm( x , f) vanish 
for t < 0. Assuming no initial surges in P(x, t) and Af(x,t) at t = to, so J e (x, to) — J m (x,to) = 0, we also have 
% e (x, t ) = Xm( x ;^o) = 0. Indeed, currents are due to the motion of charged, massive, particles and their velocity 
cannot be changed instantaneously. This property is found in linear response expressions and also, for instance, for 
the Lorentz case. Denoting dtx(t) = X'{t)i we then obtain 

Je(x,t)=/ ds X ' e (x,t-s)-E(x,s), J m (x,t) = / d SX ' m {^t-s) ■H(xl,s). (2.3) 

Jt Jt 

As a matrix, x eiin (x, t) are assumed to be symmetric (this property explicitly holds for linear response expressions if 

the unperturbed matter Hamiltonian is time-reversal invariant). 

Remarks: 

1) The initial time to can have any value, in particular to = — oo. However, in view of the Laplace-transformed 
equations, introduced later on, to — is a convenient choice. This case is often realized in practical situations. For 
instance, in describing scattering of an electromagnetic wave-packet with bounded support from a material object, 
the wave-packet is initially, as t — > — oo, well separated from the object and, in view of the hyperbolic nature of 
Maxwell's equations, the support remains bounded and, contrary to the Schrodinger case, it takes a non-zero time for 
the support to reach the object. Thus the polarization and magnetization vanish for times smaller than some finite 
to, which we set equal to 0. 

2) It is customary Q to relate M to H, rather than B, although B is the more fundamental field. Indeed, interactions 
with atoms are in terms of the microscopic vector potential A, which is related to the microscopic B-ficld. But note 
that if the particles are in vacuum, sufficiently far away from the medium, the microscopic B-field equals the 
macroscopic one and both equal H at the particle coordinates. We note further that linear response expressions 
usually relate the magnetization to the (microscopic) B-held. 

Next we introduce the Fourier decomposition 



Xe,m( x >i) = J d\ exp[-iAt]i^ e , m (x, A), 
i> e , m (x,A) = J dtexp[i\t]x' e;m (x,t) 



1 f°° 

-J dtexp[i\t]x' e , m (x,t). (2.4) 



Since %g m (x, t) are real, we have f e , m (x, — A) = i/ e ,m(x, A) in the sense that this relation holds for each component 
of these tensors. We also assume that the system is passive. This means that the electromagnetic energy 

= \f dx{E(x, tf + H(x, t) 2 }, (2.5) 
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cannot increase as a function of time. So initial population inversions in the material system are excluded. Then 
£em(t) - £em(*o) = / dsd s £ em (s) = -/ ds / dx{ J e (x, s) ■ .E(x, s) + J m (x, s) ■ Jf (x, s)} 



= — J ds du I dx{%' e (x, s — u)iE(x., u)E{s., s) + x„(x, s — u):J?(x, u)J?(x, s)} 
= -/ d£t / /*dx{Xe(x,s-u): J E(x,u) J E(x,s)+X™(x,s-u):ff(x,w)ff(x,s)} 
= — / dx[i/ e (x, A):{ / duexp[iAu].E(x, it)}{ / ds exp[— iAs]-E(x, s)} 

+ i/ m (x,A):{ / duexp[iXu]H(x,u)}{ f dsexp[-i\s]H(x, s)}] ^ 0, (2.6) 

./to « *o 

SO 

i/ e , w (x,A) >0. (2.7) 

This result also emerges in linear response expressions if the initial density operator for the material system is a 
function of its Hamiltonian and the level population decreases with increasing energy as is the case for a canonical 
distribution. Now u etTn (x, — A) = v e , m (x, A) ^ 0, leading to 

, . f ,, sin(At) f sin(At) , . , 

Xe(x,i) = / dA — - — i/ e (x,A), Xm(x,t) = / dX — - — f m (x, A), 

x' e (x,t) - / dAcos(At)i/ e (x,A), x' m (x,t) = / dA cos(At)i/ m (x, A). (2.8) 



We introduce Laplace transforms according to (Im z > 0) 

f(z) = J dtexp[izt]f(t), f(t) = ±-Jdzexp[-izt]f(z), (2.9) 

where T is a path running from — oo to +oo at some distance S > parallel to the real axis. Then, for to ^ 0, 
performing a partial integration and using Xe,m(t) = 0, t < 0, 

/>oo y* oo 

Xe,m(x, z) = / dtexp[izt\Xe, m (*,t) = dtexp{izt]xe,m(x,t) 

J t JO 

If 00 1 f°° f f 1 

= — — / dt exp[izt]Xg m (x, i) = — — / dtexp[izt] / dAf e , m (x, A) exp[— iXt] = / dXu e ^ m {x, A)— -. (2.10) 

iz Jo %z Jo J J X — z 

Causality and passivity imply that these are the most general expressions for Xe,m( x , z) [7| (the latter, being analytic in 
the upper half plane, are so-called Herglotz functions, which can always be represented in the above form). In general 
i/ e m (x, A) is made up of integrable functions, leading to absorptive systems, and of 5-function contributions. The 
second give rise to a set of dispersive Lorentz terms. In fact we can say that causal, passive, dispersive susceptibilities 
consist of a (possibly infinite) sum of dispersive Lorentz terms 0] ■ 
Then, in the isotropic, dispersive, case, x(x, z) = x(x, z)U, 

Q n fx) 2 

I/ e>ro (x, A) = ^,„(x) 2 <5(A - K, m ), Xe, TO (x, Z) = - 2 T(J ,2 - ( 2 " U ) 

Here we recall that, disregarding x-dependencies, an absorptive Lorentz contribution is given by 



X 1 {z) = -^—^ x 7 (t) =r! 2 exp[-^]^^, Ul = J u % + 'f- n,uj Oll >0. (2.12) 

z z + i-fz — ujq 2 Wi V 4 

For uel, 

f°° 1 f 1 
Xe, m (x,w)= / dtexp iwt ]Xe,m(t) = — / dA- — i/ e m (x, A), 

Im^ m (x,w) = — i/ em (x,w). (2-13) 
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In the dispersive case, 7 = 0, 

X (cu + iO) = + i*6(u> 2 ul), X (t) = 2 ^1 (2.14) 

Here the 5-function insures the validity of the Kramers- Kronig relations although it does not play a further role • 
Negative index systems are characterized by negative real permeabilities for some real frequencies. In the Lorentz 
case this can happen at a few discrete values of z. Indeed, for a single dispersive Lorentz term, setting 

l + x(z) = -l, (2.15) 

we find 



z = ±w, u = + ^, (2.16) 
but if absorption is present, only complex solutions with negative imaginary part are obtained. 

III. TIME EVOLUTION 

In this section we extend the AFF to include magnetization. Let A C M 3 be the set containing the medium, 
so v> e _ m (x,X) vanish if x ^ A as do F 2 , F 3 , F5 and Fq below. In @ E and B were used as the electromagnetic 
components. It is possible to do so in the present situation, but choosing E and H gives somewhat more symmetric 
formulae. We introduce 



F x (x,t) 


= F(x,t), 




F 2 (x,A,i) 


= U(x) /' 

J to 


ds sin{A(i 


F 3 (x,A,t) 


- U(x) /' 


ds cos{A(i 


F 4 (x,t) 


= H(x,t), 




F 5 (x,A,t) 


= U(x) /' 

J to 


ds sin{A(i 


F 6 (x,A,i) 


- i^(x) r 

J t 


ds cos{A(i 



Then 



and 



Note that 



(3.1) 



J e (x, t) = y eLW e (x, A) • y dscos{A(t — s)}Fl(x, s) = y dAi/ e (x, A) • F 6 (x, A, s), 

J m (x, *) = y rf Ai/ m (x, A) • y ds cos{A(t - s)}F 4 (x, s ) = J dXv m (x, A) • F 3 (x, A, s), (3.2) 

5 i F 1 (x,t) = 5 x xF 4 (x,t)- y dAiv e (x,A) •F 6 (x,A,i), 

t F 2 (x,A,i) = AF 6 (x,A,t), 

d t F 3 (x,A,i) = U(x)F 4 (x,t)- AF 5 (x,A,t), 

d t F 4 (x,t) = -9 X x F(x,t)-y dAiv m (x,A) ■F 3 (x,A,i), 

ftF B (x,A,t) = AF 3 (x,A,t), 

S t F 6 (x,A,i) = U(x)F x (x,*)-AF 2 (x,A,t). (3.3) 



F 2i3 , 5i6 (x,A,to) = 0. (3.4) 



In condensed notation, 

d t F(t) = -iK-F{t). 

Let (X\u e ^ m ) = i/ e>m (x, A), e be the Levi-Civita symbol and p = — id x so ( e • p) • f = id x x f). Then 



F 2 
F 3 

F 5 
\F 6 J 



, K 



/ e-p -i{v e \\ 

ooo o o a 

i\ A (x) -iX 

-e-p -i(v m \ 

iX 

ViU(x) -iX / 



K 

^ m e 



■em 

K™ o 



Thus, as in the dielectric case, K is symplcctic. Let now 

£e{t) = \ j dxj dAiy e (x,A):{F 2 (x,A,t)F 2 (x,A,t) + F 6 (x,A,t)F 6 (x,A,t)}, 
£m(t) = l / dx / dXu ™{*, X):{F 3 (x,X,t)F 3 (x,X,t)+F 5 (x,X,t)F 5 (x,X,t)}. 
Then, with £ em (t) as given by Eq. (2.5), 

dt£ em (t) = - J dx{J e (x,t) • Fi(x,t) + J m (x,t) • F 4 (x,t)}, 

d t £ e (t)= Jdxj dXu e {x, X):F 6 {x,X,t)F 1 (x,t) = / dxJ e (x,t) • Ft(x,t), 

d t £ m (t) = J dx j dXu m {x, X):F 5 (x,X,t)F i (x,t) = J dxJ m (x,t) ■ F 4 (x,t), 



so, with 



we have 



(3.5) 



(3.6) 



(3.7) 



$£(t) = o, 

i.e., £(i) is conserved in time. At this point we note that the standard expression for the conserved energy is 



(3.8) 

(3.9) 
(3.10) 



d t £ = j dx{E(x,t)-d t D(x,t) + H(x,t)-d t B(x,t)} =d t £ em (t) + J dx{E(x,t)-J e (x,t) + H(x,t)-J m (x,t)}, (3.11) 



so the two expressions agree. 
We introduce the inner product 



{F,G) = j dx[Fi(x)-Gi(x)+f4(x)-G 4 (x)]+y dx J dXv e {x, A):[F 2 (x, A)G 2 (x, A) 
+ F 6 (x,A)G 6 (x,A)] + Jjx J dXv m (x, A): [F 3 (x, A)G 3 (x, A) + F 5 (x, A) G 5 (x, A)] , 



(3.12) 



which defines the Hilbert space 



Then 



£ = ©$=i£j, K 1 =K i = L^R^dx:^), 

JC 2 =JCe = L 2 {A,dx;C 3 )® L 2 (R,v e dX), K 3 = JC 5 = L 2 {A, dx; C 3 ) ® L 2 (R, u m dX). 

|| Fit) \\ 2 = (F(t),F(t))=2£, 



(3.13) 
(3.14) 
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is conserved in time. In the Appendix it is shown (for simpler notation the isotropic case is considered) that K is 
sclfadjoint in K under some mild conditions on the susceptibilities ( \x' e m ( x , 0)| ^ c < oo), which we assume to hold 
from now on. In fact in 



K — Kq + Ki, 



(3.15) 



where 



Kn = 



V 







-e • p 





e-p \ 

iX 

o o o -a o 


iX 
-iX / 



Ki = 

















-i(v e 


\ 































i 
















-i{v m \ 

































\i 

















/ 



(3.16) 



Ki is a bounded selfadjoint operator. Thus we are dealing with a unitary time evolution on 1C, 

W(t) = exp[-tKt]. 



(3.17) 



As mentioned before it implies that the time evolution of the electromagnetic fields is properly defined. In case the 
electromagnetic fields are square integrable at the initial time, this remains true at all later times (note the passivity 
condition above). This is not evident in the original formulation which contains time convolutions. 

Remark: Note that passivity is not required to obtain the conservation of £{t). Since £ e m(0) ^ and £ e (0) = 
£m(0) = we still have £(t) = £(0) ^ 0. However, without it, £ em (t) may increase in time and £ e (t) and £ m (t) 
become negative for t > (i/ eim (x, A) may no longer be non-negative). We still can introduce the inner product 
(F, G) but the associated norm is also no longer non-negative definite. 

For dispersive systems v e . m (x, A) becomes a sum of 5- functions. Let us assume that only one dispersive Lorentz 
contribution is present in both %'s and that the medium is homogeneous and isotropic over A. Then t/ e . m (x) = 
^e,ml^i(x), so in £2,3,5,6 the x- integration is over A and 

Xe(x,i) = n 2 e \ A (x f mXet , Xei^t) = fieU(x)cOsA e i, 



x, 

sin X m t 

X m 



> Xm( x ,*) = ^„U(x)cOsA m t, 



Xe(x,z) 



A? 



rU(x), Xm(x,z) 



n 2 

L„ 



rU(x), 



(3.18) 



so 



x, s . 



J e (x,t) = figU(x) [ dscos{X e (t - s)}F(x, s) 
J t 

J m (x,i) = fi^U(x) / dscos{A m (i- s)}ff(j 
In this case, with Fl = E, F4 — H and 

F 2 (x,t) = U(x) / <issin{A e (t- s)}F(x,s), F 3 (x,t) = U(x) / rfs cos{A m (i - s)}ff(x, s), 

F 5 (x,i) = U(x) / ds sin{A m (i - s)}JT(x,s), F 6 (x,i) = U(x) / ds cos{A e (t - s)}F(x, s), 



(3.19) 



once more, 



where now 



9 t F(t) = -*K • F(t), 



(3.20) 
(3.21) 



/ 











e p 
























iX e 













«U(x) 









e p 
























2 A m 











V 


iU(x) 


-«A e 











/ 



(3.22) 
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whereas K.2,6 and £3,5 reduce to L 2 {A 1 £7 2 dx;C 3 ) and L 2 (A, f2^dx;C 3 ), respectively. In case there are more dispersive 
Lorentz terms in the susceptibilities the number of auxiliary fields increases accordingly. 

IV. LAPLACE-TRANSFORMED FIELDS 

The equations of motion can equivalently be expressed in terms of Laplace transforms. Setting t = 0, we obtain 

F(z) = i[z - K}- 1 • F(0) = iR(z) • F(0), Imz > 0. (4.1) 

From this the relations for the various components of F{z) can be obtained in terms of those of F(0) by projecting 
upon the appropriate subspace, see the Appendix. However, a direct approach involves less calculations. Thus, 

-izD(x, z) - D(x, 0) = -i(e ■ p) • ff (x, z), 

-izB(x, z) - B(x, 0) = +i(e ■ p) • £(x, z). (4.2) 

Since 

£>(x, z) = e(x, z) ■ E(x, z), B(x, z) = /x(x, z) ■ H(x, z). (4.3) 
we obtain, noting that D(x, 0) = E(x, 0) and £f(x, 0) = B(x, 0), 

L e (z) • £(x, z) - g e (x, z), L m (z) ■ ff (x, z) - g m (x, z), (4.4) 

where 

L e (z) = z 2 e(x, z) + (e • p) • /x(x, z)" 1 • (e • p), 
L m (z) = z 2 /x(x, z) + (e • p) • e(x, z)" 1 • (e • p), 
g e (x, z) = *zF(x, 0) + i(€ • p) • {/x(x, z)- 1 • if (x, 0)}, 

g m (x, z) = izH(x, 0) - i(e • p) • {e(x, z)" 1 • E(x, 0)}. (4.5) 
Here L e (z) and L m (z) are the electric and magnetic Helmholtz operators. Let now 

R e (z) = l e {z)-\ R m (z) = L m (z) _1 . (4.6) 

Then 

E(x, z) = R e (z) • g e (x, z), J?(x, z) = R m (z) • g m (x, z). (4.7) 

Note that 

<9 X • L e (z) • E(x, z) = z 2 <9 x • {e(x, z) • F(x, z)} = z 2 <9 x • £>(x, z) = izd x ■ E(x, 0) = 0, (4.8) 

as it should be. 

We can make the identification (see the Appendix) 

PiR(z)Pi = zR e (z)Pi. (4.9) 

Since the left hand side is a bounded operator, it follows that R e (z), Imz > 0, has a closed densely defined extension, 
which is in fact bounded and the same is true for R m (z). Note that, since e(x, z) = e(x, —z) and /z(x, z) = /z(x, — z), 

R e (z)* = R e (— z). Next we introduce the Green's functions 

G^"(x,y,z)-(x|R e ^(z)|y). (4.10) 

They are square integrable in x and y, respectively, are analytic in the open upper half plane and have the following 
further properties 

G(x,y,z) = G(y,x, -z) f , ^- j dzG(x, y, z) = 0, J dzzG(x, y, z) = -5(x - y)U. (4.11) 
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Now 

L e,m . G e < m ( x ,y,z) = <J(x-y)U. (4.12) 
Note that, although <9 X • L e • E(x, z) = 0, this is a special case since 

9 X • 14 • y dyG(x, y, z) • h(y) = 5 X • h(x), (4.13) 

which need not vanish for general h. 

The spatially piecewise constant situation is the case where M 3 = {UjA4.,}U{mterfaces}, with the Mj's disjoint open 
sets separated by sufficiently regular interfaces (so boundary conditions can be imposed) , whereas the susceptibilities 
are constant over Mj, 

e(x, z) = £j(z), fi(x, z) = Hj(z), xeMj. (4-14) 

V. LAYERED SYSTEMS 
A. General 

In the sequel we only consider the electric Green's function and we drop the superscript e. We also assume that 
the system is isotropic (in the anisotropic case the reduction to an expression featuring the scalar Green's functions 
for s and p polarization is not possible in general). Here we consider the situation that the A4j's are a number of 
layers parallel to the XiX 2 -plane. Then the permeabilities only depend on x 3 , 

e(x, z) = s(x 3 ,z), n(x, z) = fj,(x 3 , z). (5.1) 

We exploit the translational invariance in the X\ and ^-directions. Let k = (k\ 1 &2, k 3 ) and 
(fci,fc 2 ,0) -L e 3 , 

G(x,y,z) = (2tt)~ 2 J dncxpl-iK-ix 1 ^ - y ± )]G K (a;3, y 3 , z), 
G K (x 3 ,y 3 ,z) = (2tt)~ 2 / rfx- L exp[iK-(x" L - y ± )]G(x, y, z), 

where m can be E, B, g e '™\ etc. Then 

E(x,t) = (27r) _1 J dzexp[-izt]E(x,z) 

= {2n)~ 3 J dzexp[-izt] J dKexpl-in-tx 1 - - y- 1 )] J dy 3 G K (x 3 , y 3 , z) ■ g K {y 3 , z). (5.3) 
Omitting the subscript 3 in x 3 , etc., from now on, E K (x, z) e W = L 2 (R, dx: C 3 ), 

Sk(v, z) = izE K (y, 0) - (in + d y e 3 ) x 1 B K (y, 0), 
L K -E K (x,z) =g K (x,z), (5.4) 

and 

G K (x,y,z) = (xlL^ly), L K ■ G K (x, y, z) = S(x - y)U. (5.5) 
L K is obtained from L e by replacing p by k + pe 3 = ne K + pe 3l p = —id x . Denoting 

((x, k, z) 2 — z 2 e(x, z)fi(x, z) — k 2 , (5.6) 
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we obtain 

r ((x, n, z) 2 1 

U = {- — -, r p— ■ rp}e 3 x e K e 3 x e K 

(j,{-x.,z) 

+ {z 2 e(x, z) - p 1 p}e K e K + ^lUl^l e 3 e 3 + p K e K e 3 + —-^—pe 3 e K 

M X J Z ) M X > Z ) AH X J Z ) M X J z ) 

= L*+L£, (5.7) 

where is the s-polarization part (the term with e 3 x e K e 3 x e K ) and L£ , the p-polarization part, the remainder. 
The corresponding decomposition for G K is 

G K (x,y,z) = G s K (x,y,z) + & K {x,y,z) = (x\(L s K + (x\(L* (5.8) 

It is customary to consider the scalar Green's functions associated with the electric and magnetic s-polarization 
parts. However, the latter is transverse, whereas G£ also contains a longitudinal component. In addition, in obtaining 
atomic radiative decay rates, the full tensorial expression for the Green's function is required and there is no simple 
relation between (the transverse part of) and the magnetic s-polarized Green's function. Thus we calculated G£ 
in the Appendix with the result 

G K (x,y, z) = G s (x, y, z, k) + G p (x, y, z, k), 
G s (x, y, z, k) = G s (x, y, z, n)e 3 x e K e 3 x e K , 

G p (x, y, z, k) = (e K + -^-^d x e 3 )(e K - -^-^d y e 3 )G p (x, y, z, k), (5.9) 

where G s and G p satisfy 

{z 2 e{x, z) - p Z £<yX,Z ] p}G p {x 1 y, z,n) = 5{x-y), { K ' ^ - p 1 p}G s (x,y, z, n) = S(x - y). (5.10) 
C(x,k, z) z ^{x,z) fJ-{x,z) 

In order to obtain G p and G s we have to supplement these differential equations with the boundary conditions at 
an interface. Since d x (e/( 2 )d x G p must make sense, d x G p must exist, so we can choose G p to be continuous in x. 
In addition (e/( 2 )d x G p must be differentiable, so it must also be continuous in x. Similarly we find that G s must 
be continuous in x, as well as [i~ l d x G s and again the same applies with x and y interchanged. These boundary 
conditions can be shown to correspond to the usual boundary conditions for D and E. In addition, in view of the 
square integrability in x and y, there are no exponentially increasing contributions for layers that extend to x = ±oo. 



B. Two half spaces filled with different materials 

We consider the situation where the half spaces x > and x < are filled with media characterized according to 



We denote 



C+(k, z) 2 = z 2 e+(z)[i + {z) - k 2 , £-(«, z) 2 = z 2 S-{z)fj,-(z) - k 2 , 

K±i ^ z) = 2i^Kz)' L± ^ z) - w^Jr (512) 

and introduce the Fresnel reflection coefficients 

r P = e i + - £+C - ,r s = ^;-/ +<: - . (5.13) 

Using square integrability in x and y and the boundary conditions on the interface, we obtain in the usual way 

G p (x,y,z) = K + {exp[iC+\x - y\] - r p exp[i( + (x + y)]}9(x)9(y) + — {cxp[i( + x - iC-y]0(x)6(-y) 

iz z e+C- + 

+ exp[-iC_a; + K+y]0(-x)6(y)} + *r_{exp[<_|a; - y\] + r P exp[-*C-(a; + y)]}0(-x)O(-y), (5.14) 
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and 

G s (x,y,z) = L + {exp[i( + \x-y\]+r s exp[i( + (x + y)]}9(x)6(y) - i — + ~ {exp[iC+a; - i(-y]6(x)6(-y) 

+ exp[-i(_x + iC, + y]6(-x)6(y)} + L_{exp[i(_\x - y\] - r s exp[-iC_(ar + y)]}6(-x)6(-y), (5.15) 

from which the x and y derivatives of G p (x 7 y, z), etc., present in G p (x, y, z, k), can be obtained in explicit form. 

The case where the region x > consists of vacuum and the initial state is contained in this region is of particular 
interest. It applies to the situation where an electromagnetic wavepacket in vacuum is travelling towards the medium. 
Then, labeling vacuum quantities with the subscript and deleting the subscript — for quantities associated with the 
medium, 

K+(k,z) = K (k,z) = ^ , L+(k,z) = L (k,z) 



2iz 2 ' ' ' 2«Co(k, z)' 

Co(«, z) = r p0 = r s0 = ^-f (5.16) 



so, for y > 0, 



G p (x,y,z) = K {exp[i( \x - y\] - r p0 cxp[< (x + y)]} Q(x) + -\ ^ exp[-<x + i( y]6(-x), 



G s (a;, y, z) = L {exp[iCo|z - y\] + r s0 exp[iC (a: + y)]}0{x) - i M exp[-<x + iCot/M-a;). (5.17) 

VI. THE NIM SITUATION 

We continue our investigation of the half space case, assuming that 

e{z) = »{z) = \ (6.1) 

z — CJ n 



the dispersive Lorentz case, and rewrite 



1 = C-e(o = C-e( = C-e(o 1 1_ 

((k, z) + s(z)(o(k, z) ( 2 -e 2 ( 2 {z 2 e 2 -n 2 )-e 2 (z 2 -n 2 ) n 2 e-le + l 



Note that this expression can become infinite if e 2 = 1. In case e = 1 we are back to the vacuum case and £ — e(o = 0. 
But if £ = (j, = -1, the NIM case, 



z = uj± = ±uj, Cj = \/ + ^^ 2 ) (6-3) 



and $(z) can become infinite. Now 

1 =1 z2 -^ M Z)= C-^Co ^ 2 -^o (64) 

e(z) + l 2(z-cj)(z + cj)' lj (*-w)(* + w) 2 K 2 {e(z) - 1}' k ' 

so we encounter poles in z — ±w. Next we study the behavior of <&(z) = $(w + i<5) as S I 0. We start with ((z). In 
order to obtain the square root in the limit 5 I we must know the signs of C( w ) 2 an d cb£(cj) 2 in 

C(w + *<5) 2 = CM 2 + itfQXM 2 + 0(<5 2 ), (6.5) 

where 

n 2 \ 2 , „ / r! 2 \ / r! 2 2cj n 2 ft 2 



«*> a =*T-^J - k2 '^ 2=2 H 1 -^J1 1+ ^ + (^ifJ- (6 - 6) 

Since both quantities have definite parity it suffices to consider the case uj ^ 0. We note that C(0) 2 = ~ ft2 an d that 
C(w) 2 increases to +oo as w approaches ui . Then it decreases again to reach the value — k 2 for cj = cj = \/u!q + O 2 . 
Beyond this value it increases again to tend to +oo as cj — > +oo. Thus C( w ) 2 has three zero's, cj a £ (0, cjq), cj;, e (cj , cj) 
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and ui c > Co. Since ((u>) 2 — corresponds to a third order equation in uj 2 these are the full set of zero's. 9 w C(w) 2 
vanishes in ui = 0, then tends to +oo as uj reaches uoq, where it switches sign and increases to in Co, whereupon it 
remains positive and eventually tends to +oo. Denoting 



w 2 1 



ft 2 



U!^ — UJr 



we obtain (+ indicates that a quantity is positive, — that it is negative) 











(+UJQ, +UJb) 


(+w a , +ujq) 


(0,+w a ) 


CM 2 


+ 






+ 


+ 




Iiii((k,w + i5) 2 


+ 


+ 






+ 


+ 




p(uj) 


ip(u) 


ip(uj) 


-p(uj) 


p(uj) 


ip(uj) 


uj e 


< -w c 


(-U} c ,-U)) 


(-U), -UJb) 


(-Ul b ,-UJo) 


(-w , -uj a ) 


(-Wa,0) 


((k,uj) 2 


+ 






+ 


+ 




ImC(/t, uo + iS) 2 






+ 


+ 






C(k,uj) 


-p(uo) 


ip(u) 


ip(u) 


p(uo) 


-p(w) 


ip(u) 



In comparison, with po{uj) = \J\uo 2 — k 2 |, so p{Co) = po(Co) 



(6.7) 



(6.8) 



to e 


> K 


(0,+k) 


(-K,0) 


< -K 


Co(k,w) 




ipo(w) 


ipo(w) 


-PoM 



(6.9) 



Since uj = y to 2 + ^fl 2 we have uj € (lo ,uj). For Co > k, ((k,Cj) 2 = Co 2 — k 2 > and cD e (wo,w;,), C( k jw) = —p(Co) = 
— \/Cb 2 — k 2 , whereas for Co < k, ((k,Cj) 2 < and uj e (uJb,Co), ((k,uj) = ip{Co). Similar results follow for — Co resulting 



Co e 


> K 


(0,+k) 


({k,+Co) 


-p(Co) 


ip(Co) 


C(«,-w) 


+p{Co) 


ip(Co) 


Co(k,+w) 


+p{Co) 


ip(Co) 


Co(k,-w) 


-p(Co) 


ip(Co) 


C-(w) - £(w)Co(w) 





2ip(uj) 


C-(-w)-£(-w)Co("W) 





2ip(uj) 


CM+e(w)CoH 


-2p(Co) 





C(-w)+e(-w)Co(-w) 


+2p(Co) 






(6.10) 



For Co > k and y > we obtain 



p(Cb) 



G P {x, y, ±Co) = ±^-r^{cxp[±ip(uj)\x - y\]6(x) + cxp[±ip(Co)(x + y)]6(-x)} 
G s {x, y, ±Co) = ± ^ {exp[±ip(Co)\x - y\]9(x) + exp[±ip(Co)(x + y)]6(-x)}. 



(6.11) 



On the other hand, for Co < k, $>(z) becomes infinite for z = ±Co since C(k, z) + s(z)C,o(n, z) vanishes in ±w. Now 



$(z) 



C-eCo 



p(Co)fl 2 



{z-Co)(z + Co) 2k 2 { £ (z)-1} 



Ain 2 (z — Co)(z + Co)' 



(6.12) 
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so we encounter poles in ±Co. Thus 

G p (x,y,z) = K \exp[i( \x-y\] + - exp[iC (a; + y)]\ 6(x) + 1 ^ exp[-i(x + iCoy]6(-x) 

= K exp[i( \x - y\]6(x) + — — --{ko(( ~ e( ) exp[i( (x + y)]8(x) + -^( ( exp[-i(x + i( y]6(-x)} 

~ %T exp[-p(w)|a; - y\]6{x) + 4^2^ K _ + {<»P[-P(")(s + y)]6{x) + exp[+p(ui)(x - y)]0(-x)} 

G s (x,y,z) = L (exp[iCo|a; - y\] - j exp[iC, Q {x + y)] 1 6{x) - i M exp[-iC,x + iC,oy]6{-x) 
I s, + Mo J 4 + M4o 

= L exp[iCo|a; - v\]8(?) - — - — — {A (C - MCo) exp[iC (z + y)]0{x) + i[iexp[-i(x + iC, y]6(-x)} 

C + fKo 

Z ^ _ T7m exp[-p(u)\x - y\}0{x) + i P ^^ 7 * , ... {exp[-p(£)(s + 2/)]0(af) + cxp[+p(d))(a; - y)]0(-x)}. 

2p{Ld) AlK z [Z — 0-> ){Z + CO) 

(6.13) 

Hence, in the reflection case x, y > 0, and for Co > k, 

G p (x,y,±ui) = ±^j{exp[±ip(Co)\x - y\], G s (x,y,±Co) = ± \ {exp[±ip(Co)\x - y\], (6.14) 
Liz^ lipyio) 

where the term responsible for reflection is absent, i.e., there is no reflection at the frequencies ±w for which e(z) = 
fi(z) = -1. 

On the other hand, for Co < k, x, y > 0, 



z^±c p{Co) p(Co) 3 n 2 1 

G p {x,y,z) ~ — exp[-p(a;)|a;-y|]+ H. _ qXF z + m exp[-p{cj){x + y)}, 



G s (x,y,z) z ^~ u ___L_- e xp[-p(w)|a; - y\] + — — r r exp[-p(ui)(x + y)], (6.15) 

2p{io) 4«r (z — w)(z + w) 

so now the reflection term is still present but we encounter the damped behavior, typical for the evanescent situation. 
Next we consider refraction (transmission into the lower half space). Here y > > x, and for Co > k, 



G p (x,y,±Co) = exp[±ip(Co)(x + y)}, G s (x,y,±u) = ± . 1 exp[±ip{u) (x + y)} . (6.16) 

Alto Lipyio) 



Now 



C(w) ' 

leading to 



d x G p (x,y 7 z) ==f——, (6.17) 



G K (x,y,z) z ^ ±— — -exp[±ip(w)(y + a;)]{— ^-(e K =F -7Tre 3 )(e K =p -7Tre 3 ) + e 3 x e K e 3 x e K }. (6.18) 

iip\LO) LO z p(Ld) p{0J) 



whereas in the vacuum case 

G K (x,y, ±Co) = ± exp[±ip(Cb)(y - x)}{— r^-(e K ± — — e 3 )(e K ± -7rre 3 ) + e 3 x e K e 3 x e K }. (6.19) 

2tp{LO) LO z p(lOj p{UJ) 

Comparing the two we note that x has changed to —x and n/ p(Co) to —n/p(Co), showing the anomalous behavior found 
earlier for NIM systems. This behavior becomes more direct in a scattering formalism where it would show up in the 
corresponding scattering amplitude. However, setting up a scattering formalism, although a straightforward matter 
using the auxiliary field approach (the dielectric case was treated earlier in Q), involves a substantial amount of 
bookkeeping. This is mainly due to the existence of two different scattering channels for reflection and transmission. 
The scattering situation is somewhat less complicated for a single NIM layer, where the transmitted wave eventually 
is moving in vacuum again. 
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If Co < k, for y > > x, 

G p {x,y,z) ° exp[-p(&)(y - z)], G s (x,y,z) ^ ^ exp[-p(c2))(y - z)], (6.20) 

once more showing evanescent behavior. In retrieving _E(x,i), the pole contributions in the Green's function give 
rise to terms oscillating in time according to exp[±zd)£], so no damping occurs in the time dependence, a property 
observed earlier by Pendry [2j for the case of a single layer. 



VII. DISCUSSION 
A. Summary of results 

We started off with a system characterized by general causal, passive, susceptibilities X e (x, t) and Xm(x, t) and 
showed, using the auxiliary field approach, that E(x,t) and H(x,t) have a proper time evolution. If they are square 
integrable at the initial time this remains true at all later times so possible singularities are square integrable and the 
electromagnetic energy remains finite. We then specialized to layered systems using a Laplace transformed formalism. 
We expressed the Green's function, which is a tensor, in terms of the two scalar functions G s (x, y, z) and G p (x, y, z). A 
consequence of the auxiliary field setup is that the Helmholtz Green's function is square integrable in both coordinates. 
In particular this is true in the evanescent case. 

We then studied the special situation where one half space is vacuum and the other filled with a medium. Restricting 
ourselves to scalar permeabilities given by a single dispersive Lorentz term, 

n 2 

e(z) = n(z) = 1 - Imz > °! i 7 - 1 ) 

Z — LU 

which take on the value —1 for z = ±lj = ±(wg + ^n 2 ) 1 / 2 , we then evaluated G s (x, y, z) and G p (x, y, z) as z approaches 
these values and found, for the non-evanescent case, a typical NIM behavior, the reflected field vanishes whereas the 
transmitted field behaves anomalously, in accordance with earlier results on NIM systems. As is well known, although 
n 2 = e(±d))/i(±c2>) = 1, the NIM case is different from the vacuum situation as is seen from the refractive behavior. 
But the absence of reflection in the radiative regime is shared by both. 

On the other hand the Green's function has poles at ±a) in the evanescent situation. They lead to oscillating terms 
proportional to exp[±zwt] and the reflective part of the Green's function no longer vanishes. There is no damping 
in the temporal behavior as noted earlier by Pendry Q. This is obvious since the former only occurs in absorptive 
media, the spatial fall-off occurring in evanescent situations is completely unrelated to the temporal decay found in 
the absorptive case. 

Note that in the present setup the NIM case is a special situation occurring for two discrete frequencies ±w among 
a whole set where no NIM behavior takes place. The general expression for the permeabilities of a dispersive, non- 
absorptive, system for equal e and fi is a (possibly infinite) sum of dispersive Lorentz terms, 

£ (z) = ^) = l-E^V' ( 7 - 2 ) 

„ Z W 0n 

and in this case there is a larger set of frequencies w n for which e = fi = — 1. But between these frequencies s and 
\i vary wildly, in particular there is always a won (for which there is no transmission) between two subsequent u) n 's. 
This may spoil the idea of obtaining an interval for which this relation is approximately valid. Indeed, if in z = uj + iS, 
u approaches won, then 

e[z) ^ (7.3) 
2u! 0n d 



In the half space case this results in 



^°0, izl^ s ^-i, (7.4) 



so G K (x, y, z) vanishes in this limit if x < and y > 0. Thus the situation is opposite to the NIM case. In that case 
there is no reflection, whereas here the transmission vanishes (perfect reflector) . 
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We did not consider absorption in the NIM case. This is straightforward to do along the same lines but since the 
medium extends over a half space, the transmitted field will die out. The two poles ±u> now acquire a negative 
imaginary part, so the Green's function remains finite for all frequencies and the reflection term in the Green's 
function no longer vanishes. We intend to come back to this situation for the single layer case, where a transmitted 
field, although attenuated, is still present. Moreover this case can fairly easily be treated in terms of a scattering 
formalism, leading to scattering amplitudes for reflection and transmission, which should show NIM behavior in the 
dispersive case for appropriate frequencies. 



B. Discrete eigenvalues of K, surface modes and radiative atomic decay 

We found earlier that the Helmholtz Green's function of a NIM system had poles in z = ±lj. Usually poles in a 
Green's function originate from discrete eigenvalues of the original operator and this is precisely what happens here. 
In the more general case that A e ^ X m and fl e ^ Cl m the condition e(z)fj,(z) = 1 again gives the solutions z = ±w but 
now 

, x 2 n 2 + x 2 n 2 + n 2 n 2 

-2 /v e m 1 ' x m e 1 e m in r\ 

w = n|To|r ■ ^ 

It is straightforward to show that in the half space case ±u> are eigenvalues of K with associated eigenfunctions 
proportional to 



exp[iK ■ x^] exp[— \J n 2 — w 2 |a;3|], k > a), (7.6) 

i.e., surface modes. Since oj does not depend on k, a general theorem about direct integral decompositions tells us 
that ±w are discrete eigenvalues of K with infinite degeneracy. This can also be seen directly since by superposition 
we can construct an infinite orthonormal set of square integrable functions {f„(x)}, 

f„(x) = j dK Pn{ K)eMiK-^M-V^\**\]. (7-7) 

Here the question arises as to what happens if there are whole sets of dispersive Lorentz terms in the susceptibilities. 
Then there will be sets of solutions {±u„}, but the condition k > ui needs refinement. 

Excited atoms in vacuum decay by photon emission providing there is no selection rule forbidding the transition. 
It is well known that the radiative decay constant 7„ ac changes if we no longer have vacuum. Examples are atoms 
in a cavity or near a material such as a dielectric. This is caused by alterations in the field modes relative to the 
vacuum case. In simple situations, making the dipole and isotropic approximation, 7 can be expressed in terms of 
the imaginary part of the trace (Tr) of the Green's function as a matrix, 

7 — ImTrG(x, x,w tr ), (7.8) 

where x is the atomic position and uj tr the atomic transition frequency. In case the permeabilities are frequency- 
independent, ImTrG(x, x, ui) is the local density of states. In [6[ the above result was obtained for dielectrics. The 
same procedure, involving quantization of the fields, can be used in the present case with the same result (the case 
of an atom embedded in a magnetodielectric material was considered by the Jena group Je}). Here we make a few 
remarks about the layered case. We have 

ImTrG(x, x,wt r ) — J dK, ImTrG K (a;, x, ui tr ), (7-9) 

and consider the dispersive half space situation with the atom in vacuum close to the interface. Since the Green's 
function for x,y > is the sum of a vacuum and reflective part, we can write 7 = ^ vac + 7 re /- At the NIM 
frequency, w tr = Co, there is no reflection in the propagating regime so in the K-integral only n > Co can contribute. 
Since Im G K (x, x, ujt r + iO) becomes infinite for cotr — w, due to the presence of the 5-functions <5 (oj ± lu), we obtain 
an infinite result for 7 re /- Such an infinite local density of states was also encountered in [10(, where a so-called 
perfect corner reflector was considered. However, this result is incorrect. Upon quantization, an excited atom with 
uj tr = Co can decay radiatively but also transitions to the above bound states of K are possible. The latter process is 
an oscillatory one and we expect decay but modulated by oscillations. However, there are no infinities. 

Another point is that the situation is highly anisotropic so the isotropy approximation, which is used to obtain Eq. 
(7.8) becomes doubtful. Clearly the radiative decay problem needs further study. 
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C. Fixed frequency model 



In the introduction we mentioned that using a simple model with fixed, i.e., frequency independent e = fi = — 1 
can give rise to problems in calculating the Green's function. Indeed, our formalism indicates that, in the evanescent 
case, there are poles in ±w, so the Green's function becomes infinite and the simple model breaks down. As we have 
seen, in retrieving E(x,t) a finite result emerges and the responsible frequency integration mechanism is absent in 
the model. In addition there is a second flaw. With 

e = n = E (x 3 )=»(x 3 ) = \ +1 > X3 >°' (7.10) 

[ -1, £3 < 0, 

the conserved energy 

£ = IJ dx{e (a*) -E(x, tf + n(x 3 )B(x, tf} (7.11) 

is no longer positive definite and we can no longer base an inner product and associated Hilbert space formalism on 
this quantity. We still can use the inner product 

(f,g)= J dxf(x)-£(x). (7-12) 

but now the candidate for the generator of the time evolution, determined by the usual boundary conditions, is no 
longer selfadjoint, so the existence of a time evolution comes into question. Of course it is possible to introduce a 
Krein space with inner product based on £ but this does not solve this problem. An alternative is to restrict K to the 
eigenspaces associated with ±u), which leads to a correct time evolution. Then we have a fixed frequency model but 
now the part of K, relevant for radiative decay, is missing. 

D. Response to an external source 

We consider the time evolution due to an external source. In general the source is given by charge and current 
densities p ext (x, i) and J ext (x, t), which are related by the conservation law 

dtPextfct) +d x • Jext(x,t) = 0. (7.13) 

We assume that the source quantities vanish for t ^ to, so the same holds for the fields. Then 

d t D( x ,t) = d x x H(x,t) - J ext (x,t), d t F(t) = -iK ■ F(t) - G(t), (7.14) 

where G\(t) = J ext (x,t), whereas its other components vanish. Since F(to) = 0, Duhamel's formula gives 

F(t) = f dsexp[-iK{t-s)]-G(s). (7.15) 

We are interested in the behaviour of F(t) for large t. This depends on the nature of the spectrum of K. We assume 
that K does not have singular continuous spectrum, so 

K = ]TA„P„ + J \E ac (d\) = Y, XnPn + Yl J d\X\u Xa )(u Xa \. (7.16) 

n n a 

At this point we set t = and assume that 

G(t) = f(t)G , Go e K, (7.17) 

where the Fourier transform f(u>) in 

f(t) = J dcjcxp[-iu>t]f(u}), (7.18) 
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is a smooth function of u. Then 



P n -F(t) = f dscxp[-iX n (t-3)]f(s)P n -Go, 
Jo 

P ac -F(t) = [ rfsV / d\exp[-i\{t-s)]f(s)\u Xa )(u Xa \G ), 
Jo a J 

(u Xa \P ac -F(t)) = f dsexp[-iX(t-s)]f(s)(u Xa \G Q ), (7.19) 
Jo 



and 



P n -F(t)= [ dsexp[-iX n (t- s)]P n ■ [ du; exp[-iws]G(uj) = exp[-i\ n t] [ d^ eXp ^ A " — ^ -f(w)P n ■ G 

Jo J J l \*n — £*>) 

= cx V [-iX n t] / ^ 1 ~ 6XP[ ~ W] /(A„ + V -)P n ■ Go ' ~°° J exp[-iA„t]/(A„)P„ • G , (7.20) 
7 — iv t 2 

whereas 

(u AQ |P ac .F(t)) '-2° |exp[-iAt]/(A)<u A a|Go>, II P«c ■ f(t) H 2 '^ (|) 2 ]T / rfA|/(A)| 2 |(u Act |G )| 2 . (7.21) 

a 

We see that for large times P n -F(t) oscillates at the frequency u> n , a familiar situation for musical instruments excited 
by a pulse (playing a piano, ringing a bell). We have seen that in the dispersive half space case K has eigenvalues ±u). 
Hence E(x, t) has contributions that oscillate at these frequencies. Whether or not other contributions eventually die 
out depends on the nature of /(A) and (u Xa \Go). If the two have disjoint supports in A, this will indeed be the case. 
However, the continuous spectrum of K typically covers the whole real axis. 

Actual sources always have a finite bandwidth (although it can be quite small as for single mode laser sources). 
In practice monochromatic sources (which are usually simplified to a point source) are often considered. Thus if, for 
instance, f(t) = f smtd t 7 then /(A) = (2i)~ 1 f {S(X — uj ) — 5(u + uj )}. Now P„ • F(t) vanishes if A„ ^ ±uj and 
diverges for A„ = ±cj - Also 

{u Xa \Pac ■ t (t)) = — cxp[-i\t\{ — — }{u Xa G ), (7.22) 

2t i(\ — ujq) t(X + luq) 

and P ac ■ F{t) diverges unless (u A;Q |Go) vanishes in a neighborhood of ±o; . Such divergent behavior is typical for un- 
damped systems driven by a harmonic force. The external source situation with a monochromatic point source is some- 
times used as the starting point for the calculation of Green's functions. Then, with G(x, t) — Go exp[— iojot]5(x — xo), 
taking Fourier transforms, 

[w - K] • F(lu) = i8(uj - lu )S(x - x )G . (7.23) 

However, u — K docs not have an inverse, to being in the spectrum of K. This problem is avoided by using Laplace 
transforms, in which case 

F{z) = i[z - K]- 1 ■ G(z), Imz > 0. (7.24) 

Returning to the half space case, we conclude that for large times E(x, t) has oscillating contributions at the frequencies 
±ui but that there is also a additional contribution associated with the continuous spectrum of K. 



E. Final remarks 



The philosophy of our approach is to study special properties that occur for specific frequencies, the NIM case 
being the primary example. But, as already noted, in the half space case with a single scalar Lorentz contribution, 
transmission tends to if u> — > luq. For this to happen, e{z) and /i(z) need not be identical, as can be seen from the 
expression for the Green's function. This feature remains valid if there are more Lorentz contributions present. It 
suggests the experimental study of reflection properties as a function of frequency. Dips will occur if a NIM situation 
can happen (absorption will prevent obtaining an exact zero) and maxima at the frequencies won in Eq. (7.2). In the 
idealized single Lorentz case ujq and ui can both be obtained in this way. 



18 



Although we concentrated on the situation where the electric and magnetic susceptibilities were given by the 
same single dispersive Lorentz form, the situation where they consist of different sets of Lorentz contributions can 
also give rise to NIM situations, the fundamental requirement being the existence of some frequencies Q for which 
e(d>) = jtx(cD) = — 1. 

An important question is in how far such systems can be realized. On a microscopic level it seems not to be possible 
to obtain this property (for some further considerations on the susceptibilities of magnetodielectric systems, see fill]). 
On a larger scale, small structures, involving split rings and other configurations, lead to a more favorable situation. 
But if the electric and magnetic modes are coupled, complications arise (von Neumann's non-crossing rule). Another 
point is that, due to the employed manufacturing methods, fabricated materials are not isotropic. But this may 
change in the future. There is an extensive experimental literature concerning the fabrication of such devices A 
further complication is often the occurrence of losses. The latter can spoil the delicate effects essential for "perfect 
lenses ". In principle such lenses can consist of a single NIM slab. In the present work we did not consider this 
case, a simple form of a layered system. However, in a quite recent publication, Collin [12] made a precise analysis 
of this specific case. He took into account field contributions not considered by Pendry [2| and the effect of fields 
switched on for only a finite time interval and found that taking these into account spoils the perfect lens behavior. 
We also encountered such contributions in the half space case, the background part in addition to the pole terms in 
the Green's functions. 

Losses will blur the NIM behavior originating from pole contributions in a dispersive case, since the poles now acquire 
an imaginary part. This raises the question if adding some gain can improve the situation. Typically losses arise from a 
coupling of the electromagnetic field to material modes that have continuous spectrum. As an example, if the material 
consists of a single atom, this coupling is the one to atomic continuum states (ionization). In macroscopic media, 
among other possibilities, couplings to phonon modes and the occurrence of Forster processes can cause absorption. 
In such situations electromagnetic energy is converted to material modes where the energy leaks away in space. It 
will be clear that compensating all losses by means of adding gain (this would convert the system to a dispersive, 
non-absorptive one) will not be possible. But compensating loss at a few specific frequencies, for instance by pumping 
the system to create level inversions in the material subsystem, may be feasible. Here the gain must be controlled 
precisely in order to maintain passivity and avoid undesirable instabilities. This situation was analyzed by Stockman 
1H) 0, w ho applied a causality argument to the square of the refraction function n(z) 2 = e(z)/j(z). He concluded 
that negative refraction cannot be accomplished by adding gain. 

Finally we mention an approach based on space-time transformations leading to a simpler set of field equations 
but in a curved space-time frame TJ]. In particular the detailed work by Leonhardt and Philbin, using methods 
borrowed from general relativity, should be mentioned. So far this approach is restricted to frequency-independent 
£ and /i. Taking the auxiliary field formalism as a starting point it may be possible to extend it to the general 
frequency-dependent case. 
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Appendix A: Selfadjointness of K 



The idea is to split K into a zero order part Kq and a perturbation Ki, 



where 
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o + Ki, 
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Then Ki 
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domain V(Kq). 



19 



Proof: For notational simplicity we give the proof for scalar susceptibilities and the absorptive case (so v e ^ m (x, A) 1 / 2 
are properly defined, the dispersive situation must be handled slightly differently). Let f € 2?(K ). Then 





/gl\ 






/ 


-i J d\v e (x, A)f 6 (x, 


A) \ 




g2 















g3 




f = 




if 4 (x) 






g4 








-i J d\v m {x, A)f 3 (x 


A) 




g5 















Vg 6 y 






I 


ifi(x) 


/ 



(A3) 



Now 



gi II? 



dx J d\v e (x, A)f 6 (x, A) J dfiiy e (x, n)t e (x, 
= J dx J d\v e (x, X) 1 / 2 J d^ e (x,^ 2 iy e (x, A) 1 / 2 f 6 (x, A) • v e (x, n) 1/2 t e (x, M ) 

dx[ J d\v e {x,\t' 2 [j '^(x^)] 1 / 2 ^ ^dA Z / e (x,A)|f 6 (x,A)| 2 ] 1 /2 [ y (i/UI/e ( X)M )|f 6 ( X)/x )|2 ] i 



j dxj d\v e {x,\) I dA^ e (x,A)|f 6 (x,A)| 2 = J dx X '(x,0) J d\v e {x, A)|f 6 (x, A)| 2 <|| X '(0) 



loo || f6 \\e^ d II fb 116 



Also 



|| g 3 \\l= jdxj dA^„(x,A)|f 4 (x)| 2 < d || f 4 || 2 , 
and similar for the other components. Thus 

|| Kx f ||< Vd\\ f ||, 

so Ki is a bounded selfadjoint operator and hence K is selfadjoint with domain D(Kq).W 

Remark: Note that the proof does not require i/ 6jTO (x, A) to be non- negative or even real. However, if this is not the 
case the inner product on K. is altered and the norm no longer non-negative. 



(A4) 
(A5) 

(A6) 



Appendix B: Projections of R(z) 



We assume the susceptibilities to be scalar and consider Pi[z — K] x Pi = PiR(z)Pi, where Pj projects upon the 
j th component of f € fC, Pji = fj.We have 



[z-K]- 1 = [z + K][z 2 -K 2 }- 1 = [z + K] 



[z 2 — K em ■ K me ] 1 



K] 



[z 2 -U e ]- 1 







^h + Xe( x >°) -(^e|A -ie-p(v m \ 

-A A 2 

\-ie-p A 2 + |0)(^ m | 



, H r 



^ h o + X„(x,0) -(v m |A ie-p(i/ e | N 

-A A 2 

\*e-p A 2 + |0)(^ e |/ 



so 



P 1 [z-K}- 1 P 1 =zP 1 [z 2 -H e ]- 1 P 1 . 

According to the Feshbach projection formula with A an operator and P = 1 — Q a projector, 

A- 1 = [QAQ^Q + {P - [QAQJ^QAP^HP - PAQIQAQ]' 1 }, 
PA- X P = Q P P, Q P = [PAP - PAQ(QAQ)~ 1 QAP}- 1 , 
QA~ X Q = GqQ, Gq = [QAQ - QAP(PAP)~ 1 PAQ}~ 1 . 



(Bl) 



(B2) 



(B3) 



In our case P = Pi = 1 — Qi = 1 — P 2 — P3, A = z 2 — H e , and we obtain 

P x \z 2 - H e ]- 1 P 1 = [z 2 - h - X ' e (x,0) - PiH e Q![z 2 - QiH e Qi] _1 QiH e Pi] _1 Pi. 

Here 

P^Q^z 2 - QiHeQi^GhHePj = P!H e P 2 [z 2 - P 2 H e P 2 ]- 1 P 2 H e P 1 + P!H e P 3 [z 2 - PaHeP^PaHeP! 

= (z, e |A 2 [z 2 - A 2 ] _1 |0) - e ■ p(v m \[z 2 - A 2 - |0)(z, m |]- 1 |0) • e • p, 



and 



leading to 



so 



Similarly 



(^ e |A 2 [z 2 - A 2 ]- 1 ^) = - X ' e (x, 0) - z 2 Xe (x, z), 

(^ m \[z 2 -x 2 -\o)(^ m \r 

= -l + [l-(^ m |[z 2 -A 2 ]- 1 |0>]- 1 
= -1 + /i(x,z) _1 , 



P l [z 2 -U e ]- 1 P 1 

= [z 2 e(x, z) + (e- p)/x(x, z)" 1 • (e • p)] _1 Pi 

= Re(*)Pl, 



Pi [2 — K] _1 Pi = zR e {z)P 1 . 
Pi[z- K] _1 P 4 = zR m (z)P 4 . 



Appendix C: Decomposition of the Green's functions for the layered case 

We express R^ ,m (z) in terms of the inverses of scalar operators. Using the Feshbach formula, Eq. (B3), with 
L^ m (z) and 

P s = e 3 x e K e 3 x e K , 

Qs = U - P s = e K e K + e 3 e 3 , 



we find, noting that 



(e 3 x e K ) • (e • p) = ne 3 - pe K , 
(e • p)-(e 3 x e K ) = -ne s +pe K 
e ■ p) • Q s 

= (e • p) • (e K e K + e 3 e 3 ) 
= (e 3 x e K )(«e 3 - pe K ), 



etc., that (we skip the subscript k for brevity) 



P s -L e -P s = {^-p- P }P s , 
Ps ■ L e • Q s = 0, 

Q s ■ l e (z) -Qs = i; 

= [z 2 s(x,z)\J Q - (pe K - /ce 3 ) 1 (pe K - «e 3 )] • Q s 

/J,(x, z) 



and similar for L™ (interchange e and (i), where 

£(x, k, z) 2 = z 2 e(x, z)(j,(x, z) — k 2 . 



Thus 



where, with 



Now let 



Since 



R e > m (z) = Rf m (z) + R e p ' m (z), 



Rl(z) = [z 2 e - !^-pip}-\ R™( z ) = £ -p-p}- 

Rt(z) = Rt(z)P s , RT(z) = RT(z)P s , 

R p (z) = [z 2 eU Q - (pe R - ne 3 )-(pe K - Ke 3 )] _1 Q s 

M 

R™(z) = [z 2 fiU Q - (pe K - Ke 3 )i(pe K - Ke 3 )] _1 Q 6 



R*(z) = Ae R e K + Be K e 3 + Ce 3 e K + De 3 e 3 . 
Lp • Rp = Q s = e K e K + e 3 e 3 , 



we find by comparing coefficients that 



Let 



f K 

(z 2 e- p-p)A+p-C = 1, 
(j, [i 

C l D +-pB = l 1 

1 pA +^C = 0, 
H fj, 

1 K 

(z 2 e- p-p)B +p-D = 0. 



2 2 



Then 



R e p (z) = [z 2 e-p—p]-\ R™(z) = [z 2 ,,-p^p\ 
A = R e p (z), B = -R e p (z) P ^, C = -^pR L p (z), 

2 

f-l K> n Z £ i Ki 

D = ^ + ^Pl Z £ - P^TPl *>£2> 

with similar results for the magnetic case. Hence 

R e (z) = R e s (z)e 3 x e K e 3 x e K 

+ (e K - ^pe 3 )R e p (z){e K - p^e 3 ) + ^e 3 e 3 , 



R m (z) = RT(z)e 3 x e K e 3 x e K 
+ (e« - ^pe 3 )R™{z)(e K - p-^e 3 ) + ^e 3 e 3 . 
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Note that 

w = lSs RT{z)p ^ (C14) 

so Rp(^) can be expressed in terms of the scalar magnetic s-polarized Green's function but the vectors in front and 
behind are quite different. From the above expressions the corresponding Green's functions, introduced in Sect. V, 
now follow. 
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